%% Imperfect competition
clc;
clear;

N_M=4;

N_I=10;
N_U=100;

lambda_I =1;
lambda_U =1;

T=10;
tau_V=1;



tau_eps=2/T*ones(T,1);
tau_eps(1)=0.5;


tau_eta=2*T*ones(T,1);
tau_eta(1)=10;



theta_bar=1;
V_bar=1;

X_0=0;


Theta_I=theta_bar;
Theta_U=theta_bar-Theta_I;

theta_I_bar=Theta_I/N_I;
theta_U_bar=Theta_U/N_U;




o_I_V=zeros(T,1);
for j=1:T
    o_I_V(j)=1/(tau_V+sum(tau_eps(1:j)));
end
K_I_V=o_I_V.*tau_eps;

sigma_I=zeros(T,1);     
sigma_I(1)=1/tau_V+1/tau_eps(1); 
for j=2:T
    sigma_I(j)=o_I_V(j-1)+1/tau_eps(j);  
end

Sigma_I=zeros(2,2,T);
for t=1:T
    Sigma_I(1,1,t)=sigma_I(t);
    Sigma_I(2,2,t)=1/tau_eta(t);
end



% mu=zeros(T,1);
% mu(T-1)=lambda_I*o_I_V(T);
% for t=T:-1:2
%     mu(t-1)=mu(t)+lambda_I*sigma_I(t)*(K_I_V(t))^2;
% end

mu= lambda_I*o_I_V;

h=lambda_I./tau_eps;

Sigma_U=zeros(T,1);
Sigma_U(1)=1/tau_V+1/tau_eps(1)+h(1)^2/tau_eta(1);

o_U_V=zeros(T,1);
o_U_X=zeros(T,1);
o_U_VX=zeros(T,1);

o_U_V(1)=(tau_eps(1)^-1+h(1)^2*tau_eta(1)^-1)/(tau_V*Sigma_U(1));
o_U_X(1)=1/tau_eta(1)-h(1)^2/(tau_eta(1)^2*Sigma_U(1));
o_U_VX(1)=(1/tau_V)*(h(1)/tau_eta(1))/Sigma_U(1);

K_U_V=zeros(T,1);
K_U_X=zeros(T,1);

K_U_V(1)=tau_V^-1/Sigma_U(1);
K_U_X(1)=-h(1)*tau_eta(1)^-1/Sigma_U(1);

for j=2:T
    Sigma_U(j)=o_U_V(j-1)+1/tau_eps(j)+h(j)^2/tau_eta(j);
    
    o_U_V(j)=(1/tau_eps(j)+h(j)^2/tau_eta(j))*o_U_V(j-1)/Sigma_U(j);
    o_U_X(j)=o_U_X(j-1)+1/tau_eta(j)-(o_U_VX(j-1)-h(j)/tau_eta(j))^2/Sigma_U(j);
    o_U_VX(j)=o_U_VX(j-1)-(o_U_VX(j-1)-h(j)/tau_eta(j))*o_U_V(j-1)/Sigma_U(j);
    
    K_U_V(j)=o_U_V(j-1)/Sigma_U(j);
    K_U_X(j)=(o_U_VX(j-1)-h(j)/tau_eta(j))/Sigma_U(j);
end

H_I=zeros(4,4,T);
for t=2:T
    H_I(1,1,t)=1;
    H_I(2,2,t)=1;
    H_I(3,3,t)=1-K_U_X(t)*mu(t-1);
    H_I(3,2,t)=K_U_X(t)*mu(t-1);
    H_I(4,4,t)=1;
end

F_I=zeros(4,2,T);
for t=1:T
    F_I(1,1,t)=K_I_V(t);
    F_I(2,2,t)=1;
    F_I(3,1,t)=K_U_X(t);
    F_I(3,2,t)=-K_U_X(t)*h(t);
end

F_U=zeros(3,1,T);
for t=1:T
    F_U(1,1,t)=K_U_V(t);
    F_U(2,1,t)=K_U_X(t);
end








% the values of the parameters for the last period

gamma_I=zeros(T,1);
gamma_I(T)=lambda_I*o_I_V(T);

gamma_U=zeros(T,1);
gamma_U(T)=lambda_U*o_U_V(T);

w_I=zeros(T,1);
w_I(T)=(N_M/(N_M+1))*(gamma_I(T)/N_I)/(gamma_I(T)/N_I+gamma_U(T)/N_U);

w_U=zeros(T,1);
w_U(T)=(N_M/(N_M+1))*(gamma_U(T)/N_U)/(gamma_I(T)/N_I+gamma_U(T)/N_U);

delta_I=zeros(T,1);
delta_I(T)=1;

delta_U=zeros(T,1);
delta_U(T)=1;

g_I=zeros(T,1);

g_U=zeros(T,1);

f_I=zeros(T,1);

f_U=zeros(T,1);

rho_I=zeros(T,1);
rho_I(T)=1;

rho_U=zeros(T,1);
rho_U(T)=1;

m_I=zeros(T,1);
m_I(T)=gamma_I(T);

m_U=zeros(T,1);
m_U(T)=gamma_U(T);




C_I=zeros(4,1,T);
C_I(3,1,T)=-w_I(T)*mu(T);

C_U=zeros(3,1,T);
C_U(2,1,T)=w_U(T)*mu(T);

M_I=zeros(4,4,T);
M_I(1,2,T)=1;
M_I(2,1,T)=1;
M_I(2,2,T)=-mu(T);
M_I(3,3,T)=(w_I(T)*mu(T))^2/gamma_I(T);

M_U=zeros(3,3,T);
M_U(2,2,T)=(w_U(T)*mu(T))^2/gamma_U(T);



Xi_I=zeros(2,2,T);
Xi_I(:,:,T)=eye(2)/( eye(2)/Sigma_I(:,:,T)+ lambda_I*transpose(F_I(:,:,T))*M_I(:,:,T)*F_I(:,:,T));

Xi_U=zeros(T,1);
Xi_U(T)= 1/(1/Sigma_U(T)+ lambda_U*transpose(F_U(:,:,T))*M_U(:,:,T)*F_U(:,:,T));


H_I_P=zeros(4,1,T);
H_I_P(1,1,T)=(1-w_I(T))*delta_I(T)+w_I(T)*delta_U(T);
H_I_P(2,1,T)=-mu(T)*((1-w_I(T))*delta_I(T)+w_I(T)*delta_U(T));
H_I_P(3,1,T)=(1-w_I(T))*g_I(T)+w_I(T)*(g_U(T)+delta_U(T)*mu(T));
H_I_P(4,1,T)=(1-w_I(T))*f_I(T)+w_I(T)*f_U(T);

H_U_P=zeros(3,1,T);
H_U_P(1,1,T)= w_U(T)*delta_I(T)+ (1-w_U(T))*delta_U(T);
H_U_P(2,1,T)= w_U(T)*(g_I(T)-delta_I(T)*mu(T))+ (1-w_U(T))*g_U(T);
H_U_P(3,1,T)= w_U(T)*f_I(T)+ (1-w_U(T))*f_U(T);

H_I_theta=zeros(4,1,T);
H_I_theta(1,1,T)= w_I(T)*(delta_I(T)-delta_U(T))/gamma_I(T);
H_I_theta(2,1,T)= -w_I(T)*(delta_I(T)-delta_U(T))*mu(T)/gamma_I(T);
H_I_theta(3,1,T)= -w_I(T)*(mu(T)*delta_U(T)+g_U(T)-g_I(T))/gamma_I(T);
H_I_theta(4,1,T)= - w_I(T)*(f_U(T)-f_I(T))/gamma_I(T);

H_U_theta=zeros(3,1,T);
H_U_theta(1,1,T)=-w_U(T)*(delta_I(T)-delta_U(T))/gamma_U(T);
H_U_theta(2,1,T)= w_U(T)*(delta_I(T)-delta_U(T))*mu(T)/gamma_U(T)+...
    w_U(T)*(mu(T)*delta_U(T)+g_U(T)-g_I(T))/gamma_U(T);
H_U_theta(3,1,T)=w_U(T)*(f_U(T)-f_I(T))/gamma_U(T);


H_I_theta_hat=zeros(4,1,T);

H_U_theta_hat=zeros(3,1,T);








test=zeros(T,1);
test(T)=mu(T);

SOD_I=zeros(T,1);
SOD_I(T)=gamma_I(T);

SOD_U=zeros(T,1);
SOD_U(T)=gamma_U(T);

% compute the values of the parameters recursively

for t=T:-1:2
    
    SOD_U(t-1)=lambda_U*Xi_U(t)*(transpose(H_U_P(:,:,t))*F_U(:,:,t))^2;
    
    gamma_U(t-1)=gamma_U(t)/(N_M+1)+ lambda_U * Xi_U(t)* transpose(H_U_P(:,:,t))*F_U(:,:,t)*...
        transpose(F_U(:,:,t))*(H_U_P(:,:,t)+ C_U(:,:,t)/(N_M+1));
    
    Phi_U= transpose(H_U_P(:,:,t))- lambda_U*Xi_U(t)* transpose(H_U_P(:,:,t))*F_U(:,:,t)*...
        transpose(F_U(:,:,t))*M_U(:,:,t);
    
    delta_U(t-1)=Phi_U(1,1);
    
    g_U(t-1)=Phi_U(1,2);
    
    f_U(t-1)=Phi_U(1,3)- (gamma_U(t)+lambda_U * Xi_U(t)* transpose(H_U_P(:,:,t))*F_U(:,:,t)*...
        transpose(F_U(:,:,t))*C_U(:,:,t) )*w_I(t)/N_U;
    
    
    SOD_I(t-1)=lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*H_I_P(:,:,t);
    
    gamma_I(t-1)=gamma_I(t)/(N_M+1)+ lambda_I* transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*(H_I_P(:,:,t)+C_I(:,:,t)/(N_M+1));
    
    Phi_I= transpose(H_I_P(:,:,t))*H_I(:,:,t)- lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*...
        Xi_I(:,:,t)*transpose(F_I(:,:,t))*M_I(:,:,t)*H_I(:,:,t);
    
    delta_I(t-1)=Phi_I(1,1);
    
    g_I(t-1)=Phi_I(1,3);
      
    f_I(t-1)=Phi_I(1,4)-(gamma_I(t)+ lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*C_I(:,:,t) )*w_U(t)/N_I;
    
    test(t-1)=-Phi_I(1,2)/delta_I(t-1);   % mu(t-1)
    
    
    w_U(t-1)=(N_M/(N_M+1))*(gamma_U(t-1)/N_U)/(gamma_I(t-1)/N_I+gamma_U(t-1)/N_U);
    w_I(t-1)=(N_M/(N_M+1))*(gamma_I(t-1)/N_I)/(gamma_I(t-1)/N_I+gamma_U(t-1)/N_U);
    
    
    H_I_P(1,1,t-1)=(1-w_I(t-1))*delta_I(t-1)+w_I(t-1)*delta_U(t-1);
    H_I_P(2,1,t-1)=-mu(t-1)*((1-w_I(t-1))*delta_I(t-1)+w_I(t-1)*delta_U(t-1));
    H_I_P(3,1,t-1)=(1-w_I(t-1))*g_I(t-1)+w_I(t-1)*(g_U(t-1)+delta_U(t-1)*mu(t-1));
    H_I_P(4,1,t-1)=(1-w_I(t-1))*f_I(t-1)+w_I(t-1)*f_U(t-1);
    
    H_U_P(1,1,t-1)= w_U(t-1)*delta_I(t-1)+ (1-w_U(t-1))*delta_U(t-1);
    H_U_P(2,1,t-1)= w_U(t-1)*(g_I(t-1)-delta_I(t-1)*mu(t-1))+ (1-w_U(t-1))*g_U(t-1);
    H_U_P(3,1,t-1)= w_U(t-1)*f_I(t-1)+ (1-w_U(t-1))*f_U(t-1);
    
    H_I_theta(1,1,t-1)= w_I(t-1)*(delta_I(t-1)-delta_U(t-1))/gamma_I(t-1);
    H_I_theta(2,1,t-1)= -w_I(t-1)*(delta_I(t-1)-delta_U(t-1))*mu(t-1)/gamma_I(t-1);
    H_I_theta(3,1,t-1)= -w_I(t-1)*(mu(t-1)*delta_U(t-1)+g_U(t-1)-g_I(t-1))/gamma_I(t-1);
    H_I_theta(4,1,t-1)= - w_I(t-1)*(f_U(t-1)-f_I(t-1))/gamma_I(t-1);
    
    H_U_theta(1,1,t-1)=-w_U(t-1)*(delta_I(t-1)-delta_U(t-1))/gamma_U(t-1);
    H_U_theta(2,1,t-1)= w_U(t-1)*(delta_I(t-1)-delta_U(t-1))*mu(t-1)/gamma_U(t-1)+...
       w_U(t-1)*(mu(t-1)*delta_U(t-1)+g_U(t-1)-g_I(t-1))/gamma_U(t-1);
    H_U_theta(3,1,t-1)=w_U(t-1)*(f_U(t-1)-f_I(t-1))/gamma_U(t-1);
    
    H_I_theta_hat(:,:,t-1)=H_I_theta(:,:,t-1)/(N_M+1);
    H_I_theta_hat(4,1,t-1)=H_I_theta_hat(4,1,t-1)+w_U(t)/N_I;
    
    H_U_theta_hat(:,:,t-1)=H_U_theta(:,:,t-1)/(N_M+1);
    H_U_theta_hat(3,1,t-1)=H_U_theta_hat(3,1,t-1)+w_I(t)/N_U;
    
    
    
    m_U(t-1)=m_U(t)/(N_M+1)^2+ lambda_U * Xi_U(t)* (transpose(H_U_P(:,:,t))*F_U(:,:,t))^2-...
        lambda_U* Xi_U(t)*(transpose(C_U(:,:,t))*F_U(:,:,t))^2/(N_M+1)^2;
      
    
    m_I(t-1)=m_I(t)/(N_M+1)^2+ lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*H_I_P(:,:,t)- lambda_I*transpose(C_I(:,:,t))*...
        F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*C_I(:,:,t)/(N_M+1)^2;
    
    
    C_I(:,:,t-1)= transpose(H_I(:,:,t))*C_I(:,:,t)/(N_M+1)+ m_I(t)*H_I_theta_hat(:,:,t-1)/(N_M+1)+...
        lambda_I* H_I_theta(:,:,t-1)*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*H_I_P(:,:,t)- (lambda_I/(N_M+1))*transpose( C_I(:,:,t)*...
        transpose(H_I_theta_hat(:,:,t-1))+ M_I(:,:,t)*H_I(:,:,t) )*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*C_I(:,:,t);
    
    C_U(:,:,t-1)= C_U(:,:,t)/(N_M+1) + m_U(t)*H_U_theta_hat(:,:,t-1)/(N_M+1)+...
        lambda_U* Xi_U(t)*(transpose(H_U_P(:,:,t))*F_U(:,:,t))^2*H_U_theta(:,:,t-1)-...
        (lambda_U*Xi_U(t)/(N_M+1))*transpose( C_U(:,:,t)*transpose(H_U_theta_hat(:,:,t-1))+...
        M_U(:,:,t))*F_U(:,:,t)*transpose(F_U(:,:,t))*C_U(:,:,t);
    
    M_I(:,:,t-1)=transpose(H_I(:,:,t))*M_I(:,:,t)*H_I(:,:,t)+m_I(t)*H_I_theta_hat(:,:,t-1)*...
        transpose(H_I_theta_hat(:,:,t-1))+ H_I_theta_hat(:,:,t-1)*transpose(C_I(:,:,t))*...
        H_I(:,:,t)+ transpose(H_I(:,:,t))*C_I(:,:,t)*transpose(H_I_theta_hat(:,:,t-1))+...
        lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*...
        H_I_P(:,:,t)*H_I_theta(:,:,t-1)*transpose(H_I_theta(:,:,t-1))-...
        lambda_I*transpose( C_I(:,:,t)*transpose(H_I_theta_hat(:,:,t-1))+ M_I(:,:,t)*H_I(:,:,t) )*...
        F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*( C_I(:,:,t)*...
        transpose(H_I_theta_hat(:,:,t-1))+ M_I(:,:,t)*H_I(:,:,t) );
    
    M_U(:,:,t-1)=M_U(:,:,t) +m_U(t)*H_U_theta_hat(:,:,t-1)*transpose(H_U_theta_hat(:,:,t-1))+...
        H_U_theta_hat(:,:,t-1)*transpose(C_U(:,:,t))+C_U(:,:,t)*transpose(H_U_theta_hat(:,:,t-1))+...
        lambda_U*Xi_U(t)*(transpose(H_U_P(:,:,t))*F_U(:,:,t))^2*H_U_theta(:,:,t-1)*...
        transpose(H_U_theta(:,:,t-1))-lambda_U*Xi_U(t)*transpose(C_U(:,:,t)*...
        transpose(H_U_theta_hat(:,:,t-1))+M_U(:,:,t))*F_U(:,:,t)*transpose(F_U(:,:,t))*...
        ( C_U(:,:,t)*transpose(H_U_theta_hat(:,:,t-1))+M_U(:,:,t));
    
    Xi_I(:,:,t-1)= eye(2)/(eye(2)/Sigma_I(:,:,t-1)+...
        lambda_I*transpose(F_I(:,:,t-1))*M_I(:,:,t-1)*F_I(:,:,t-1));
    
    Xi_U(t-1)= 1/(Sigma_U(t-1)^-1+ lambda_U*transpose(F_U(:,:,t-1))*M_U(:,:,t-1)*F_U(:,:,t-1));
    
    rho_I(t-1)=rho_I(t)*sqrt(det(Xi_I(:,:,t))/det(Sigma_I(:,:,t)) );
    
    rho_U(t-1)=rho_U(t)*sqrt(Xi_U(t)/Sigma_U(t) );   

end

SOC_I=find(SOD_I<0)
SOC_U=find(SOD_U<0)















% Trading Vol and spread



tau=(1./tau_eps+h.^2./tau_eta).^-1;

d_V=zeros(T,1);
d_eps=zeros(T,T);
d_eta=zeros(T,T);



for t=1:T

    
    d_V(t)= (w_U(t)/(gamma_U(t)/N_U))* ( (delta_U(t)+(g_U(t)-g_I(t))/mu(t))*o_U_V(t)*sum(tau(1:t))-...
        (delta_I(t)+(g_U(t)-g_I(t))/mu(t))*o_I_V(t)*sum(tau_eps(1:t)) );
    
    for n=1:t
        d_eps(t,n)= (w_U(t)/(gamma_U(t)/N_U))* ( (delta_U(t)+(g_U(t)-g_I(t))/mu(t))*o_U_V(t)*tau(n)-...
            (delta_I(t)+(g_U(t)-g_I(t))/mu(t))*o_I_V(t)*tau_eps(n) ) ;
        
        d_eta(t,n)= (w_U(t)/(gamma_U(t)/N_U))*( (delta_I(t)+(g_U(t)-g_I(t))/mu(t))*mu(t)-...
            (delta_U(t)+(g_U(t)-g_I(t))/mu(t))*o_U_V(t)*tau(n)*h(n) );
    end
end



        
n_M=zeros(T,1);
n_M(1)=N_M/(N_M+1);
for j=2:T
    n_M(j)=n_M(j-1)/(N_M+1);
end


D_V=zeros(T,1);
D_eps=zeros(T,T);
D_eta=zeros(T,T);

for t=1:T
    D_V(t)=-d_V(t)+ sum(n_M(t-1:-1:1).*d_V(1:t-1));
    
    for n=1:t
        D_eps(t,n)=-d_eps(t,n)+ sum(n_M(t-n:-1:1).*d_eps(n:t-1,n));
        
        D_eta(t,n)=-d_eta(t,n)+ sum(n_M(t-n:-1:1).*d_eta(n:t-1,n));
    end
end


            

sigma_Delta_sq=zeros(T,1);
for t=1:T
    sigma_Delta_sq(t)= D_V(t)^2/tau_V + sum(D_eps(t,:).^2./transpose(tau_eps))+...
        sum(D_eta(t,:).^2./transpose(tau_eta));
end
sigma_Delta=sqrt(sigma_Delta_sq);




mu_Delta=zeros(T,1);

mu_Delta(1)=Theta_U*N_M/(N_M+1) + ( w_U(1)/(gamma_U(1)/N_U)*(f_U(1)-f_I(1))-w_I(1) )*theta_bar;
for t=2:T

      mu_Delta(t)=Theta_U*N_M/(N_M+1)^t +( (w_U(t)/(gamma_U(t)/N_U))*(f_U(t)-f_I(t))-w_I(t) )*theta_bar...
          -theta_bar*sum(n_M(t-1:-1:1).*( (w_U(1:t-1)./(gamma_U(1:t-1)/N_U)).*(f_U(1:t-1)-f_I(1:t-1))-w_I(1:t-1) ) );

end



Vol=   sqrt(2/pi)*sigma_Delta.*exp( -0.5*(mu_Delta./sigma_Delta).^2 ) +...
    mu_Delta.*( 1-2* normcdf(-mu_Delta./sigma_Delta) )  ;

Vol_tot=sum(Vol);

Spread= Vol.*(gamma_I/N_I+gamma_U/N_U)/N_M;









figure(1)
plot(0:1/(T-1):1, Spread, '--');
xlabel('$t/T$','interpreter','latex')
ylabel('Expected Bid-Ask Spread','interpreter','latex')
hold on
legend('$N_M=1,T=10$', '$N_M=1,T=30$', '$N_M=4,T=10$','interpreter','latex')











